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The MetropoUs Monte Carlo algorithm with the Finite Element method applied to compute elec- 
trostatic interaction energy between charge densities is described in this work. By using the Finite 
Element method to integrate numerically the Poisson's equation, it is shown that the computing 
time to obtain the acceptance probability of an elementary trial move does not, in principle, depend 
on the number of charged particles present in the system. 
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Long ranged electrostatic interactions are particularly 
difficult to take into account correctly in computer simu- 
lation of condensed matter. At present and to this aim, 
the most adequate methods, in accuracy and efficiency, 
are the Ewald summations [1-3]. In classical electrody- 
namics, the long ranged feature of electrostatic inter- 
actions is included in the structure of the local partial 
differential equation for the electrostatic field (Poisson, 
Laplace, etc.) and the particular analytical form of these 
interactions is governed by the boundary conditions [4]. 
In Ewald sums, one uses (partial) periodic boundary con- 
ditions [5, 6] and a macroscopic boundary condition (sur- 
face term) to obtain univocally energy of systems [5-10] ; 
then the energy computed likewise is used in Metropolis 
algorithm to sample the phase space of the system [11] in 
Monte Carlo simulations or, in Molecular dynamics, one 
computes forces in similar ways. 

The Finite Element method allows to obtain numeri- 
cal approximations of partial differential equations with 
boundary conditions [12-14] ; applied to Poisson's equa- 
tion, they define Poisson solvers that may be used to com- 
pute electrostatic fields and energy. Such Poisson solvers 
have been used in ab initio electronic structure calcu- 
lations in solids [15, 16]. The Finite Element method 
combined with Ewald sums has also been proposed re- 
cently for the electrostatic interactions in slablike geom- 
etry [17]. 

In this letter, we show how to use the Finite Element 
method to obtain acceptance probability of trial moves in 
Metropolis Monte Carlo algorithms in systems with elec- 
trostatic interactions. The work presented here aims to 
describe theoretically the algorithm ; to be considered as 
an efficient alternative to Ewald methods, detailed com- 
parisons arc still to be made. The main result of the 
present work is the derivation of the electrostatic energy 
of a trial configuration from the energy of the old con- 
figuration and the variations of few shape functions that 
gives locally an estimate of the electrostatic field (Eq.l9). 
For a charge density distribution p{x) on a domain 57 of 
the space, the electrostatic potential is given by the Pois- 



son's equation as 

VV = -Anp (1) 

On boundaries P of 17 some conditions must be fulfilled 
by the potential to obtain an unique solution. These 
boundary conditions can be chosen by specificying the 
value of the electrostatic potential on P (Dirichlet prob- 
lem) or the value of the electric field (Neumann problem) 
[4]. One can also show that Eq.(l) has an unique solution 
for mixed boundary conditions : taking Dirichlet bound- 
ary conditions on P/j and Neumann boundary conditions 
on P^v with P = Pij U Pat. For mixed boundaries condi- 
tions, the electrostatic potential is constrained by 

WxeTo, V{x) = Vo{x) (2) 

and 

dV 

WxeTN, -g^h{x) = -Eoix) (3) 

where Vq and Eq arc known functions and they are called 
boundary conditions on P. Eqs. (1-3) define the so-called 
strong form of a boundary-value problem in the finite 
element method [12]. The Finite Element method starts 
by deriving a weak formulation of the boundary value 
problem ; strong and weak formulations of the boundary- 
value problem are fully equivalent. The weak formulation 
for Poisson's equation reads as follows : for any weighting 
functions (or variations) w G S, find the function V & T 
such as 

[ Vw.VV dn = An [ wp dn~ [ w n.Eo dT (4) 

where the set T of trial functions and the set S of weight- 
ing functions are defined respectively by 

r T = {veH^n) : yx eTD,v{x) ^Va{x)} 

{ (5) 

[e ^ {w e H\n) : yx e Td, w{x) = 0} 
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with H^{il), the Sobolev space of degree one [12, 14]. 
The boundary conditions given by Eq.(2) are called es- 
sential boundary conditions, they are explicitly involved 
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in the definition of T, while the boundary conditions 
given by Eq.(3) are cahed natural boundary conditions 
and they are incoporated into Eq. (4) by using the Green 
first identity in the derivation of the weak formulation. 
If the essential boundary conditions have finite values on 
each point on F^i, then the weighting functions w are 
required to satisfy homogeneous counterparts of the es- 
sential boundary conditions [12], as defined by the set S, 
Eq.(5). 

Having the weak form of a boundary- value problem, the 
Finite Element method then provides an approximation 
for the function V. The first step is to construct finite- 
dimensional representations T'' and 2'' as approxima- 
tions for T and S ; we require that T'' C T and S'' C S. 
Both sets T'' and S'* are associated to a mesh that di- 
vides n into domain elements whose characteristic length 
or volume is represented by h. In the following, we re- 
strict the discussion to Galcrkin's approximation scheme 
for the Finite Element Method of C°-class (for instance 
see Chap.2 of ref.[13] and refs.[12, 14]). 
On element domains, nodal points and shape functions, 
associated to each nodal point, are defined and numbered 
(locally and globally). We use the following notations : 
the position of a nodal point numbered i is ; the set 
of all nodal points is ui ; the set of all nodal points lying 
on F^) is called the set of nodal points not lying 

onTu is called w and the function is the shape func- 
tion associated to the nodal point numbered i. We note 
also I w I, the number of nodal points in lo. Wc require 
(j)i{xj) = Sij, where Sij is the Kroneckcr symbol. The 
nodal points are particular locations in fl where the ap- 
proximate solution given by the Finite Element method 
is the most accurate. 

For G E^, the method assumes that it can be written 
as 

w;'^(a;) =^c,0,(a;) (6) 

where c,; is a constant, i.e. the value of it;'' at node i. For 
g g^^^ j g have Ci ~ 0. The approximate 

solution of the weak form admits the representation 

V'' = v'' + Vq^ (7) 

where G S'* and result in satisfaction of the bound- 
ary condition given by Eq.(2), these functions are ex- 
panded according to 

where the vt are the unknowns values of at each node 
i and 

Vo^ix) = M^'dM^) (9) 

With this representation, the integral equation Eq.(4) 
becomes an algebraic equation with | lo \ unknowns. One 



may easily obtain 

^a((?!),,0j>j =47r(0,,p) -(0,, n.£;o)r„ 

- ^ a(0j,0fc)Vb(a;fc) 

(10) 

where we have set 

= a{4>i,4>j) = / V(t)i.V4>j dVl 
Jn 

< (0i,p) = / (j)^pdn (11) 

Jn 

((^j,n.£;o)r„ = / (tiifi.EodV 

Eq.(lO) is the matrix form of the weak formulation of the 
boundary-value problem in the Galerkin approximation 
scheme. An estimation of V{x) is therefore obtained as 
in Eq.(7) by solving the linear algebraic equations given 
by Eq.(lO) to obtain the constants Vj. With convenient 
choices of nodes and shape functions, the matrix in the 
left handed side of Eq.(lO) is sparse. 
When the boundary conditions used are periodic, for in- 
stance in ah initio solid-state electronic structure compu- 
tations [15, 16] or in Molecular simulations, for all x on 
the boundaries (F) of fi, the electrostatic potential must 
fulfill 

r V{x) = V{x + L^) 

{ (12) 
[n{x).'VV{x) = -n{x + La).'VV{x + La) 

Where the domain J7 is chosen as a parallelepipedic cell 
(orthorhombic or triclinic), h the outward normal unit 
vector and La are the lattice vectors associated with the 
periodic boundaries conditions. These boundary con- 
ditions are easily fulfilled with a clever global number- 
ing of nodal points (cf. ref.[15, 16]). Partial periodic 
boundary conditions are used to represent quasi-one and 
quasi-two dimensional systems (nanotubes, surfaces, in- 
terfaces, biological membranes, etc.) ; for these systems, 
Eqs.(12) are applied for directions having spatial peri- 
odicities, while Eqs.(2,3) for other directions. Thus, the 
method seems particularly well adapted to the represen- 
tation of systems with partial periodic boundary condi- 
tions. 

For a set of point charges {qa}a=i,N located at 
{xa}a=i.N, the charge density is given by 

N 

p{x) ^YqaS{x - Xa) (13) 

then, the electrostatic potential is given by Eqs.(7,8) with 
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Vj the solutions of Eq.(lO). More precisely, 

N 

Vji{Xa}) ^ 47r^qa^{A^^)jtMXa) 

OL—1 iGuj 

-^(A-i)j4(</),,n.£;o)r„ + J2 ai<p^,(|)k)Vo{xk)] 

(14) 

where A ^ stands for the inverse matrix of matrix A de- 
fined in Eq.(ll). The approximation for the electric field 
is E'^ = —'VV^ and, for natural boundary conditions 
(Neumann) , the energy of the configuration is given by 

W = — f \E''\'^ dx 

^ (15) 

An elementary trial move for the particle numbered a is 
Xa ^ x'g^ ; from this elementary move the trial config- 
uration of the set of charges is {xa}'. Then, the trial 
configuration is accepted or rejected according to the 
Metropolis algorithm. From Eq.(lO) and (11), the elec- 
trostatic potential for the trial configuration is given by 
the electrostatic potential of the old configuration as 

w,({a;„}') =w»({a;„})-f 47rga^(A-i)yA0j(a) (16) 
where 

Ac^,{a) = c^,{x',)~cf>,{xa) (17) 

More precisely, the approximate solution of the electro- 
static potential for the trial configuration is thus given 
by 

V''i{Xa,}',x) ^V'\{Xa,},x) 

+ATTqa M^) X1(^~^)»J^'?^J (") 

(18) 

and, from Eqs.(15) and (16), we obtain 
W =W + qaY.v,{{x^})A(b^{a) 

(19) 

This last equation is the main result of this letter. From 
Eq.(19), we may accept or reject the trial move by us- 
ing the standard Metropolis scheme with the probability 
min(l,exp(-/3(VF' - Vr))). 

If the trial move is accepted, then the new electrostatic 
potential is obtained by updating Vi{{xa}) using Eq.(16) 
; thus an another trial move may be done from the 



new configuration. If the trial move changes the loca- 
tion of more than one charge, as it is the case for poly- 
atomic molecules such as water or even more complicated 
molecules as DNA, proteins, RNA, etc., then, since equa- 
tions are linear, Eq.(16) and (19) must include a summa- 
tion over charges that arc involved in the trial move. 
In principle, such algorithm is very efficient since each 
charge is located in only one element domain and the 
shape function is different from only in domain el- 
ements that contain the nodal point Xi associated with 
the shape function 0^. If there are n nodal points in each 
domain element, then there are at most 2n values A0j (a) 
defined by Eq.(17) that are non zero and that contribute 
to Eq.(19). Once all lS.(j)j(a) are computed, to evaluate 
W' from W the number of additions and products needed 
scale as A-n? . 

For concrete implementations of this method in a stan- 
dard NVT Monte Carlo simulations, one may proceed as 
follows : 

(a) The domain 57 is defined as the simulation box, 
it is divided into domain elements ; locations of 
nodal points in each elements and shape functions 
for each nodal point are defined. Boundary condi- 
tions Eqs.(2,3) and Eq.(12) are applied. (This step 
defines the matrix A.) 

(b) The inverse of matrix A is computed and stored. 

(c) An initial configuration for charges (or charge den- 
sity) is build. (After this step, the charge density 
is defined and with step (a) all contributions in the 
right handed side of Eq.(lO) are defined.) 

(d) The algebraic equation (10) is solved by using the 
matrix A~^ . (This step gives an approximation of 
the electrostatic potential as Eq.(7).) 

(e) Electrostatic energy of initial configuration is com- 
puted by using Eq.(15). (The cpu-time of this step 
scales as N"^ where N is the number of charges in 
the domain VL.) 

(f ) Trial moves are performed and they are accepted or 
rejected according to the Metropolis algorithm by 
using Eq. (19) ; averages are accumulated. (Accord- 
ing Eq.(19), only matrix A~^ , computed in step (b), 
is needed to compute the energy of the trial config- 
uration from the energy of the old configuration ; 
thus, the cpu-time needed for each trial moves does 
not depend on N .) 

This method may also be applied to Monte Carlo simula- 
tions that sample the shape of the simulation box (NPT, 
etc.). In these implementations, steps (a) to (e) have to 
be done for each trial move of the box shape. 
In its principles, the Monte Carlo Finite Element algo- 
rithm presented in this letter has several benefits in com- 
parison to others algorithms for electrostatic long ranged 
interactions as Ewald methods, etc. 
First, and as the main benefit, the algorithm is local. The 
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summations over the pairs of particles is done only one 
time in the whole computation, at step (e). As a conse- 
quence of Eq.(19), the cpu-time needed to compute the 
acceptance probability of a trial move does not depend 
on the number of particles in the system ; this makes the 
method far more efficient than Ewald methods. 
By comparison with mesh based Ewald methods (PME, 
SPME, P^M, etc.) [2, 3], it is not necessary to assign 
charge to mesh points : each charge keeps its own iden- 
tity. 

By comparison with auxiliary field local algorithm devel- 
oped in refs.[18], one does not have to introduce a gauge 
invariant constrained auxiliary field in the Hamiltonian 
of the system to compute energy. 

Finally, another important benefit of this algorithm is 
practical. Finite Element Method is one of the most used 
method to solve numerically partial differential equa- 
tions, therefore there are a lot of very efficient codes that 
already exists. 

Nevertheless, the algorithm presented in the present let- 
ter can not yet be considered as an alternative to exist- 
ing methods ; to be so, the efficiency and accuracy of the 
Monte Carlo Finite Element algorithm has to be com- 
pared to the efficiency and accuracy of Ewald methods. 
Such a comparison is still to be made, preferencially on 
a simple system. 

There are mainly two points that strongly affect the ac- 
curacy and efficiency of the algorithm. First, the number 
of domain elements and the number of nodal points are 
critical for the convergence rate of approximate solutions 
V'^ to electrostatic potential V. These two numbers de- 
fine the size of matrix A and A~-^ : | w | x | w |. If the 
value I a; I is too large to obtain an accuracy similar to 
the one achieved in Ewald methods, then the efficiency 
may drastically fall because of processors memory access 
problems. Second, it is necessary to compute energy of 



the system as in Eq.(15) and not as 

W=^ [ p{x)V^{x)dx (20) 

Otherwise, contributions to energy from self interactions 
are lost (cf. Chap.l of ref.[4]) and comparisons with 
Ewald methods become difficult to achieve. Clearly, this 
implies that Finite Elements of C^-class have to be used 
to ensure the continuity of the electric field across bound- 
aries of domain elements ; to apply the method to molec- 
ular dynamic simulations this requirement is essential 
too. 

Equation (8), that defines the unknowns Vi, corresponds 
to element of C°-class : the approximation is continu- 
ous on r2, but derivatives of are not continuous across 
boundaries of domain elements. For elements of C^-class, 
the equation that defines the unknowns is written as 

v\x) = Y^[nM^)+v'rU^r\^)) (21) 

with the same notations as in equation (8) and (/),["^ 
the shape functions for partial derivatives (they verify 
: (l)f\xj) = and da(j)'°'\xj) = 5ij) and w-"^ the 
unknowns for derivatives. In one hand, because of equa- 
tion (21), the number of unknowns by domain elements 
of C^-class is increased by comparison to elements of 
C°-class. thus is the size of matrix A and On 
the other hand, the convergence rate of Finite Element 
Method of C^-class is extremely better than C°-class. 
Therefore some optimal choices for shape functions, 
structure of domain elements, etc. are to be made. 
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